package ru.vsu.amm.calculus;


public class Cholesky {

    
    private double precision;
    private double[] solution;
    
    private Cholesky() {}
    
    public static Cholesky process(double[][] lenta, double[] f){
        double[][] C = new double[lenta[0].length][];
        double[][] B = new double[lenta[0].length][];       
        
        for(int i = 0; i < lenta[0].length; i++){
            int minus = lenta[0].length - i < lenta.length ? lenta.length - lenta[0].length + i : 0;
            B[i] = new double[lenta.length-minus];
            C[i] = new double[lenta.length-minus];
            for(int j = 0; j < lenta.length-minus; j++){
                double Summ = 0;
                int kBoundary = Math.min(lenta.length-j-1, i);
                for (int k = 1; k <= kBoundary; k++){
                    Summ+=B[i-k][j+k]*C[i-k][k];
                }
                B[i][j] = lenta[j][i]-Summ;                
                C[i][j] = B[i][j]/B[i][0];
            }
        }
        
        // make f1 for x = (1, 1, 1, ..., 1).
        double[] f1 = f1builder(lenta);
        
        // BCx=f, By=f -> y - ?
        double[] y = new double[C.length];
        double[] y1 = new double[C.length]; //to calculate precision
        for (int i = 0; i<C.length; i++){
            double Summ = 0;
            double Summ1 = 0;
            int j = 1;
            int jBoundary = Math.min(i, C[0].length-1);
            for (; j <= jBoundary; j++){
                Summ+=B[i-j][j]*y[i-j];
                Summ1+=B[i-j][j]*y1[i-j];
            }
            y[i] = (f[i] - Summ) / B[i][0];
            y1[i] = (f1[i] - Summ1) / B[i][0];
        }       
        
        // Cx=y, x - ?
        // y like x to save memory
        for (int i = C.length-1; i>=0; i--){
            double Summ = 0;
            double Summ1 = 0;            
            for (int j = 1; j < C[i].length; j++){
                Summ+=C[i][j]*y[i + j];
                Summ1+=C[i][j]*y1[i + j];
            }
            y[i] = (y[i] - Summ) / C[i][0];
            y1[i] = (y1[i] - Summ1) / C[i][0];
        }
        
        Cholesky ch = new Cholesky();
        ch.solution = y;
        
        double abs;
        for (int i=0; i<f1.length; i++){
            if ((abs = Math.abs(1-y1[i])) > ch.precision) {
                ch.precision = abs;
            }
        }
        
       return ch;        
    }
    
    public double getPrecision() {
        return precision;
    }
    
    public double[] getSolution() {
        return solution;
    }

    private static double[] f1builder(double[][] lenta){
        double summ;        
        double[] f1 = new double[lenta[0].length]; 
        for (int i=0; i<lenta[0].length; i++){
            int minus = lenta[0].length - i < lenta.length ? lenta.length - lenta[0].length + i : 0;
            summ=0;
            for (int j=0; j<lenta.length - minus; j++){
                summ+=lenta[j][i];
            }
            int jBoundary = Math.min(i, lenta.length - 1);
            for (int j=1; j<=jBoundary; j++){
                summ+=lenta[j][i-j];
            }
            f1[i]=summ;
        }
        
        return f1;
    }
}
